usingPlotsS(a,b,θ) =a*b*sin(π/6)+b^2*sin(θ) *cos(θ)+a^2*sin(π/6- θ)*cos(π/6- θ)# アニメーションの作成@giffor b in1:0.005:2plot(x ->S(1, b, x), 0, π/6, label =false, title ="S(1, $(round(b,digits=1)), θ)", legend =:topright)end every 1
[ Info: Saved animation to /Users/shimizudan/Documents/Teaching_materials/2024/20250327tokyo-u/tmp.gif
usingPrintfusingOptim# 関数 S(a, b, θ)S(a, b, θ) =a*b*sin(π/6) + b^2*sin(θ)*cos(θ) + a^2*sin(π/6- θ)*cos(π/6- θ)# θ ∈ [0, π/6] の範囲で最大値を求める関数functionmax_S(a, b) result = Optim.optimize(θ ->-S(a, b, θ), 0.0, π/6) # 最大化なのでマイナスを最小化 θ_max = Optim.minimizer(result) S_max =S(a, b, θ_max)return θ_max, S_maxend# 各 a, b の組に対して計算a =1b_list = [1+.1*i for i=0:10]for b in b_list θ, Sval =max_S(a, b)@printf"a=1, b=%.4f → θ=%.4f, max S=%.6f\n" b θ Svalend
a=1, b=1.0000 → θ=0.2618, max S=1.000000
a=1, b=1.1000 → θ=0.3434, max S=1.109933
a=1, b=1.2000 → θ=0.4132, max S=1.239062
a=1, b=1.3000 → θ=0.4708, max S=1.385884
a=1, b=1.4000 → θ=0.5177, max S=1.548764
a=1, b=1.5000 → θ=0.5236, max S=1.724279
a=1, b=1.6000 → θ=0.5236, max S=1.908513
a=1, b=1.7000 → θ=0.5236, max S=2.101407
a=1, b=1.8000 → θ=0.5236, max S=2.302961
a=1, b=1.9000 → θ=0.5236, max S=2.513176
a=1, b=2.0000 → θ=0.5236, max S=2.732051
b_list = [1.4+.01*i for i=0:10]for b in b_list θ, Sval =max_S(a, b)@printf"a=1, b=%.4f → θ=%.4f, max S=%.6f\n" b θ Svalend
a=1, b=1.4000 → θ=0.5177, max S=1.548764
a=1, b=1.4100 → θ=0.5219, max S=1.565878
a=1, b=1.4200 → θ=0.5236, max S=1.583127
a=1, b=1.4300 → θ=0.5236, max S=1.600468
a=1, b=1.4400 → θ=0.5236, max S=1.617895
a=1, b=1.4500 → θ=0.5236, max S=1.635409
a=1, b=1.4600 → θ=0.5236, max S=1.653010
a=1, b=1.4700 → θ=0.5236, max S=1.670697
a=1, b=1.4800 → θ=0.5236, max S=1.688471
a=1, b=1.4900 → θ=0.5236, max S=1.706331
a=1, b=1.5000 → θ=0.5236, max S=1.724279
第4問
第4問・問題
Julia言語
素数パッケージ Primes.jl を利用
usingPrimes,Printff(a::Int, x::Int) = x^2+ x - afunctionN(a) k =0for n =1:a val =f(a, n)if val ≥0&&isqrt(val)^2== val k +=1endendreturn kendprintln(" a | N(a) | 4a+1 | isprime?")println("-----|------|------|----------")for a =1:120 count =N(a) val =4a +1 is_p = count ==1 ? string(isprime(val)) :""@printf("%4d | %4d | %4d | %s\n", a, count, val, is_p)end
usingPlots# 曲線C上の点を生成(z such that |z - 1/2| = 1/2)functiongenerate_C(N=500) θ =range(0, 2π, length=N)return [1/2+1/2*cis(t) for t in θ if abs(1/2+1/2*cis(t)) >1e-6] # 原点を除くendplot(1./generate_C(),aspectratio=true,xlim=(-2,2),ylim=(-2,2),label="1/z",legend=false)